Dynamic critical behavior of the classical anisotropic BCC Heisenberg antiferromagnet 
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Using a recently implemented integration method [Krech et. al.] based on an iterative second- 
order Suzuki- Trotter decomposition scheme, we have performed spin dynamics simulations to study 
the critical dynamics of the BCC Heisenberg antiferromagnet with uniaxial anisotropy. This tech- 
nique allowed us to probe the narrow asymptotic critical region of the model and estimate the 
dynamic critical exponent z — 2.25 ± 0.08. Comparisons with competing theories and experimental 
results are presented. 



I. INTRODUCTION 

Despite their relative simplicity, classical Heisenberg 
models can be used to describe the dynamic critical be- 
havior of magnetic materials whose ions have large ef- 
fective spin values. For example, RbMnFs is well mod- 
eled by an isotropic Heisenberg antiferromagnet, whereas 
MnF2 and FeF2 are good physical realizations of the 
Heisenberg antiferromagnet with weak and strong uni- 
axial anisotropics, respectively. 

The effect of uniaxial anisotropy on the critical dy- 
namics has been the subject of many theoretical 
and experimental studies 0, . Although both indicate 
that the dynamic structure factor has a dominant purely 
diffusive longitudinal component and a suppressed prop- 
agative transverse component, there is still controversy 
about the dynamic critical exponent z. While a dynamic 
scaling argument Q and mode-coupling theory Q indi- 
cate that z = 2, a combination of renormalization group 
theory and e-expansion ^ predicts that z = 2 + a/ v. Us- 
ing the Ising value v = 0.6289(8) from a high resolution 
Monte Carlo study P| and a — 0.110(5) from field theory 
0, the latter prediction is z — 2.175(10). Experiments 
on FeF2 Q have not had enough resolution to distin- 
guish between these two predictions, and on MnF2 
the exponent obtained was that of the isotropic antifer- 
romagnet, z — 1.5. Presumably this very low value is due 
to crossover effects resulting from the weak anisotropy in 
MnF2. 

Our approach to this problem is to use spin dynam- 
ics simulations with an efficient integration scheme. The 
dynamic critical exponent is estimated using a dynamic 
finite-size scaling theory and compared with the predic- 
tions by the different theories and with experiments. 



II. MODEL AND METHODS 

The Hamiltonian of the model is given by 

J2 Sr -Sr' -i?^(5/)2 (1) 
<r,r'> r 

where Sr = {S^ , S^, S^) is a three-dimensional classical 
spin of unit length, < r, r' > denotes nearest-neighbor 



pairs of spins coupled by the antiferromagnetic exchange 
interaction J > 0, and D is the uniaxial single-site 
anisotropy. We use D = 0.0591J and the inverse critical 
temperature 1/Tc = 0.478fc_B/J appropriate for MnF2 
We consider body-centered-cubic lattices with lin- 
ear sizes 12 < L < 60, corresponding to 2L^ sites, and 
periodic boundary conditions. 

The time evolution of the spins is governed by cou- 
pled equations of motion, and the dynamic structure 
factor 5(q, w), for momentum transfer q and frequency 
Lu, observable in neutron scattering experiments, is the 
Fourier transform of the space-displaced, time-displaced 
spin- spin correlation function C*''(r — r' ,t), defined as 

C^-(r-r',t) ^ (^r'W^r'(O)) - (5^(i))(5^,(0)). 

with k = x,y, or z. Because of the periodic boundary 
conditions, we can only access a set of discrete values of 
momentum transfer given by q — 2Tmq/L, where riq = 
±l,±2,...,±L/2. 

Equilibrium configurations at Tc were generated using 
a hybrid Monte Carlo method where each hybrid step 
consisted of two Metropolis sweeps through the lattice 
and eight Wolff cluster spin flips. Typically 3000 hybrid 
Monte Carlo steps were used to generate each equilibrium 
configuration, which was then used as an initial spin con- 
figuration in the integration of the coupled equations of 
motion. We use 2040 initial configurations to compute 
the averages. For the largest lattice used {L = 60), there 
are 432,000 coupled equations to integrate. These nu- 
merical integrations were performed to a maximum time 
of tmax = 400/ J, using an iterative second-order Suzuki- 
Trotter decomposition method implemented by Krech et 
al 0, with two iterations and a time step of St = 0.04/ J. 
This method consists of performing explicit rotations of 
a spin about its local effective field. The lattice is divided 
into two sublattices A and B, and a formal solution of 
the equations of motion yields y{t + St) = e'^+^)''*y(i), 
where y = {yA, ys) denotes collectively all the spins, and 
the matrices A and B are the infinitesimal generators of 
rotation of the spin configurations yA on sublattice A at 
fixed ys and of the spin configurations ys on S at fixed 
yA, respectively. The spin configurations on sublattices A 
and B are updated in succession, using the second-order 
Suzuki- Trotter decomposition of the exponential opera- 
tor El, given by e^^+B)** = e^^'^^e^^^e^^'/^ + 0{St^). A 
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similar scheme using a fourth-order decomposition of the 
exponential operator has also been tested in Ref . ^ . The 
iterative procedure is necessary because of the single-site 
anisotropyQ. 

The finite integration time introduces oscillations when 
the Fourier transform of C'^(r — r',t) is taken to pro- 
duce S'(q, Lu) and we smooth them out by convoluting 
the spin-spin correlation function with a Gaussian reso- 
lution function in frequency, with characteristic width 6uj ■ 
The dynamic structure factor thus obtained is denoted as 

The dynamic critical exponent z can be deter- 
mined from dynamic finite-size relations given by 
cc>5i(q,w)/xi(q) = Giu;L^,qL,S^L^) and 

Urn^L-'n{qL,6^L'), (2) 

where xi(l) = X!!^ '^'^(q: t^)c?'^/27r is the momentum 
dependent susceptibility and tD™ is a characteristic fre- 
quency, given by /"^^ 5|;(q, w)dw/27r = xi(q)/2. To es- 
timate z we choose the width of the resolution function to 
be = 0.015(60/L)^ so that the function n{qL,S^L'') 
in Eq. ||2J) is a constant if qL is fixed, yielding 

UJrn ~ L'""- (3) 

Because depends on z, this exponent had to be deter- 
mined iteratively. The coefficient 0.015 of 6^^ was chosen 
empirically as a compromise between effectively reducing 
the oscillations in 5(q, u) and not excessively broadening 
its structure. 



III. RESULTS AND DISCUSSION 

Our simulations, performed at the critical tempera- 
ture for MnF2 Q, show that while the transverse com- 
ponent S'^{c[,t) of the space- Fourier-transformed corre- 
lation function S'(q, t) has a short relaxation time, the 
longitudinal component (q, t) decays extremely slowly 
[see Fig.(^], and this requires that the equations of mo- 
tion be integrated to very long times. Although our tmax 
is still not large enough for the longitudinal component to 
decay to a value that is close to zero, it is a significant im- 
provement over the maximum time we could reach with 
more standard integration methods, such as the fourth- 
order predictor-corrector method, with our current com- 
puting resources. Whenever not shown, error bars are 
smaller than the symbol sizes in the figures. 

The longitudinal component of the dynamic structure 
factor, 5'^(q, w), shown in Fig.||2K), has a purely dissi- 
pative behavior, as indicated by the single central peak. 
In contrast, the transverse component S'^{q,uj) contains 
a propagative mode, indicated by the spin wave peak in 
Fig. (lib); with a possible small central peak as well. A 
comparison between Figs.(|2i) and ^p) shows that the 
longitudinal component has a much larger intensity. 
The dispersion curve i^{q) of the spin waves along the 
[100] direction is shown in Fig.JSJl, where the solid line 




100 200 300 

tJ 



FIG. 1: Longitudinal and transverse components of the space- 
Fourier-transformed correlation function. A few characteristic 
error bars are shown. 
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FIG. 2: (a) Longitudinal and (b) transverse components of 
the dynamic structure factor, with Sui = 0.015. 

represents a fit to the function lu = ujq + cq^ . The non- 
zero value of the frequency in the limit as (7 —> is charac- 
teristic of an anisotropic model. These results for the dy- 
namic structure factor are in qualitative agreement with 
theory and experimental results. 

An estimate for the dynamic critical exponent z was 
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FIG. 3: Dispersion relation for the spin wave excitation in 
^"^(q,^;) for small values of q. 




FIG. 4: Estimate of the dynamic critical exponent for differ- 
ent Hq. Analysis done with L = 30, 36, 48 and 60. 



obtained iteratively using Eq.Q for different values of 



Uq. Such estimates are denoted as Zq and are shown 
in Fig. 10} . For the weakly anisotropic system considered 
here, the onset of the asymptotic critical region occurs at 
very low values of q, accessible only with very large lattice 
sizes (and not accessible to experiments yet). Therefore, 
a better estimate for z is given by zq, obtained by ex- 
trapolating Zq to the limit q ^ 0. We fitted Zq with the 
function Zq = zq + aUq + hnq"^, using Uq = 1, 2, Uq'^^^, 
where zq, a, b, and c are fitting parameters. To check 
the robustness of the extrapolated zq, we have also per- 
formed fittings with fixed c = 2. The systematic change 
in the value of zq for both types of fittings was studied as 
the number of points included in the fittings increased, 
from Uq"^^'^ = 4 and 5, for the fittings with c = 2 and 
variable c, respectively, to n^™'*'^ = 15. For Uq"^'^^ up to 
7, the of the fittings is reasonable, and there is no clear 
systematic change in the value of zq; hence, the average 
Zq from these fittings yields z — 2.25 ± 0.08, which is the 
first numerical estimate of the dynamic critical exponent 
for this model. 



IV. SUMMARY AND CONCLUSIONS 

Because of difficulties for experiments to probe the crit- 
ical region, experimental data have not yet been able to 
distinguish between competing theories. While limited 
by finite lattice size and finite integration time, simula- 
tions offer the hope of shedding light on the differences 
between theories and experiment. Although not yet con- 
clusive, our estimate of z = 2.25 ± 0.08 is slightly larger 
than, but consistent with, the prediction by the renor- 
malization group theory. It is not consistent with mode- 
coupling theory and the dynamic scaling prediction. 

This work was partially supported by NSF grant DMR- 
0094422. Simulations were performed on the Cray T90 
and IBM SP (Blue Horizon) at SDSC. 
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